#include<poisson_construct.h>
#include<tool/libcell_kd_node_iterator.h>
#include<tool/libcell_octree_iterator.h>
#include<Mesh/libcell_iterator.h>

#include<clapack/include/libcell_clapack.h>
#include<libcell_macros.h>
double poisson_oc_node_width(OC_Node* octree)
{
    for(auto oit=octree_begin(octree);oit.it!=NULL;oit++)
    {
        if(oit.it->value!=NULL)
        {

            return oit.it->loc_max[2]-oit.it->loc_min[2];
        } 
    }

}
static inline void poisson_oc_prop_init(Poisson_OC_Prop* pop)
{
    memset(pop->normal,0,sizeof(double)*3);
    memset(pop->coeffs,0,sizeof(double)*3);
}
static inline void poisson_oc_prop_free(Poisson_OC_Prop* pop)
{
    if(pop==NULL)
    {
        return ;
    }
    free(pop);
}

void poisson_oc_node_init_from_node(OC_Node*ocn,Node* node)
{
    oc_node_init(ocn);
    if(node==NULL)
    {
        return ;
    }
    Node*nit=node;
    template_v* v=(template_v*)(nit->value);
    ocn->loc_min[0]=v->point[0];  ocn->loc_min[1]=v->point[1]; ocn->loc_min[2]=v->point[2];
    ocn->loc_max[0]=v->point[0];  ocn->loc_max[1]=v->point[1]; ocn->loc_max[2]=v->point[2];
    for(;nit!=NULL;nit=(Node*)(nit->Next))
    {
        v=(template_v*)(nit->value);
        for(int i=0;i<3;i++)
        {
            v->point[i]>ocn->loc_max[i]?(ocn->loc_max[i]=v->point[i]):(v->point[i]<ocn->loc_min[i]?(ocn->loc_min[i]=v->point[i]):(i=i));
        }
    }
    double max=ocn->loc_max[0]-ocn->loc_min[0];
    max=((ocn->loc_max[1]-ocn->loc_min[1])>max?(ocn->loc_max[1]-ocn->loc_min[1]):(max));
    max=((ocn->loc_max[2]-ocn->loc_min[2])>max?(ocn->loc_max[2]-ocn->loc_min[2]):(max));
    ocn->loc_max[0]=ocn->loc_min[0]+max;ocn->loc_max[1]=ocn->loc_min[1]+max;ocn->loc_max[2]=ocn->loc_min[2]+max;
   
    ocn->value=node_copy(node);
}

double fx_dot_fx[3]={11/20.0,13/60.0,1/120.0},fx_dot_Dfx[3]={0, 5/12.0,1/24.0}, fx_dot_DDfx[3]={-1,1/3.0,1/6.0};

//fx_dot_DDfxi=DDfxi_dot_fx
//fx_dot_Dfxi=-Dfxi_dot_fx

double* compute_value_according_basis_function(double*x,OC_Node*octree,int depth)
{
    double * re=(double*)malloc(sizeof(double)*3);
    memset(re,0,sizeof(double)*3);
    double sum=0;
    for(auto oit=octree_begin(octree);oit!=octree_end(octree);oit++)
    {
        Poisson_OC_Prop* pop=(Poisson_OC_Prop*)(oit.it->prop);
        if(pop!=NULL&&oit.it->depth==depth)
        {
            double value=oc_poisson_FF(x,oit.it);
            if(value>0.0000001)
            {
                for(int j=0;j<3;j++)
                {
                    re[j]+=value*pop->normal[j];
                }
                printf("value:%lf\n",value);
            }
            sum+=value; 
        }
    }
    printf("summmmm:%lf\n",sum);
    return re;
}
Node*  compute_support_functions_nodes(double * x ,KD_Node*kdtree,double w )
{
    Node* re=NULL;
    Vertex *v=(Vertex*)malloc(sizeof(Vertex));
    Vertex_init_(v); 
    v->point_size=3;
    v->point=(double*)malloc(sizeof(double)*3);
    v->point[0]=x[0];v->point[1]=x[1];v->point[2]=x[2];    
    re=kd_tree_find_nearest_sphere(v,3*w,kdtree);

    free_Vertex(v);
    return re;
}
double * compute_value_according_basis_function1(double*x,KD_Node* kdtree,double w)
{
    double* re=(double*)malloc(sizeof(double)*3);
    memset(re,0,sizeof(double)*3);
    double sum=0;
    Node* n=compute_support_functions_nodes(x,kdtree,w);
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        OC_Node* ocnode=(OC_Node*)(v->prop); 
        Poisson_OC_Prop* pop=(Poisson_OC_Prop*)(ocnode->prop);
        double value=oc_poisson_FF(x,ocnode);
        for(int j=0;j<3;j++)
        {
            re[j]+=value*pop->normal[j];
        }  
        sum+=value;
    }
    printf("summmmm:%lf\n",sum);
    return re;
}


void points_cloud_compute_normal_using_clapack(Node* n,KD_Node*kdtree,int num)
{
    //printf("libo\n");
    int i=0;
    double *normal=NULL;
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
    
        template_v* v=(template_v*)(nit->value);
        Node* n1=kd_tree_find_n_adjacent_verticesn(v,kdtree,num);
        double* axis=compute_pca_result_using_clapack(n1);  
        normal=(double*)malloc(sizeof(double)*3);
        normal[0]=axis[2];normal[1]=axis[5];normal[2]=axis[8];
        v->prop=normal; 
        free(axis);
        free_node(n1);
        i++;
    } 
    //printf("libo\n");
 
} 

void compute_octree_normal(OC_Node* octree,int depth)
{
    for(auto oit=octree_begin(octree);oit!=octree_end(octree);oit++)
    {
        if(oit.it->value==NULL)
        {
            continue;
        }   
        Poisson_OC_Prop* pop=(Poisson_OC_Prop*)malloc(sizeof(Poisson_OC_Prop));
        poisson_oc_prop_init(pop);
        for(Node* nit=oit.it->value;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_v* v=(template_v*)(nit->value);
            double* normal1=(double*)(v->prop);
            pop->normal[0]+=normal1[0]; pop->normal[1]+=normal1[1]; pop->normal[2]+=normal1[2]; 
        }
        normalize(pop->normal,3);   
        oit.it->prop=pop;
    }   

    //以平均值的方法计算八叉树空节点的法向。
    for(auto oit=octree_begin(octree);oit!=octree_end(octree);oit++)
    {
        if(oit.it->depth==depth&&oit.it->value==NULL)
        {
            Poisson_OC_Prop* pop=(Poisson_OC_Prop*)malloc(sizeof(Poisson_OC_Prop));
            poisson_oc_prop_init(pop);  
            for(int i=0;i<8;i++)
            {
                OC_Node* ocnode=oit.it->parent->children[i];
                Poisson_OC_Prop* pop1=(Poisson_OC_Prop*)(ocnode->prop);
                if(pop1!=NULL)
                {
                    pop->normal[0]+=pop1->normal[0];pop->normal[1]+=pop1->normal[1];pop->normal[2]+=pop1->normal[2];
                }
            }   
            normalize(pop->normal,3);           
            oit.it->prop=pop;
        }
    }
}


Mesh* poisson_transfer_octree2vertices(OC_Node*octree,int depth)
{
    Mesh* re=(Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);
    double p[3]={0}; 
    // double max=0;
    for(auto oit=octree_begin(octree);oit!=octree_end(octree);oit++)
    {
        if(oit.it->prop==NULL)
        {
            continue;
        }
        //printf("depth:%d\n",oit.it->depth); 
        for(int j=0;j<3;j++)
        {
            p[j]=(oit.it->loc_min[j]+oit.it->loc_max[j])/2.0;
        }
       // max=oit.it->loc_max[0]-oit.it->loc_min[0];
        template_v* v= re->create_vertexv(re,p,3);
        v->prop=oit.it;
    }
    //printf("max :%lf\n",max);
    return re;
}

double poisson_compute_b_i(template_v* v,KD_Node*kdtree,double w)
{
    double re=0;     
    Node* n=kd_tree_find_nearest_sphere(v,4*w,kdtree);
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v1=(template_v*)(nit->value);
        OC_Node* ocnode=(OC_Node*)(v1->prop); 
        Poisson_OC_Prop* pop=(Poisson_OC_Prop*)(ocnode->prop);
        int weizhi[3]={0};
        for(int j=0;j<3;j++)
        {
            double value1=fabs(v1->point[j]-v->point[j]);
            weizhi[j]=(value1<0.5*w?(0):( value1<1.5*w?(1):(value1<2.5*w?(2):(-1))));
        }
        if(weizhi[0]==-1||weizhi[1]==-1||weizhi[2]==-1)
        {
            continue;
        }
        for(int j=0;j<3;j++)
        {
            double value=pop->normal[j];
            for(int k=0;k<3;k++)
            {
                if(k==j)
                {
                    value*=(fx_dot_Dfx[weizhi[k]]/w);
                } 
                else
                {
                    value*=fx_dot_fx[weizhi[k]];
                }
            }
            re+=value;
        } 
    }
    // printf("summmmm:%lf\n",sum);
    return re;
}

double poisson_compute_A_i_j(template_v* vi,template_v* vj ,double w)
{
    double re=0;
    int weizhi[3]={0};
    for(int j=0;j<3;j++)
    {
        double value1=fabs(vj->point[j]-vi->point[j]);
        weizhi[j]=(value1<0.5*w?(0):( value1<1.5*w?(1):(value1<2.5*w?(2):(-1))));
    }
    if(weizhi[0]==-1||weizhi[1]==-1||weizhi[2]==-1)
    {
            return re;
    } 
    for(int j=0;j<3;j++)
    {
        double value=1;
        for(int k=0;k<3;k++)
        {
            if(k==j)
            {
                value*=(fx_dot_DDfx[weizhi[k]]/(w*w));
            } 
            else
            {
                value*=fx_dot_fx[weizhi[k]];
            } 
        }
        re+=value;
    } 
    return re;
}






